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ABSTRACT 

We report the detection of eight bright X-ray bursts from the 6.5 s magnetar IE 1048.1—5937, during a 2013 July 
observation campaign with the Nuclear Spectroscopic Telescope Array. We study the morphological and spectral 
properties of these bursts and their evolution with time. The bursts resulted in count rate increases by orders of 
magnitude, sometimes limited by the detector dead time, and showed blackbody spectra with kT ~ 6-8 keV in 
the 790 duration of 1-4 s, similar to earlier bursts detected from the source. We find that the spectra during the tail 
of the bursts can be modeled with an absorbed blackbody with temperature decreasing with flux. The burst flux 
decays followed a power law of index 0.8-0. 9. In the burst tail spectra, we detect a ~13 keV emission feature, 
similar to those reported in previous bursts from this source as well as from other magnetars observed with the 
Rossi X-ray Timing Explorer. We explore possible origins of the spectral feature such as proton cyclotron emission, 
which implies a magnetic field strength of B ~ 2 x 10 15 G in the emission region. However, the consistency of the 
energy of the feature in different objects requires further explanation. 
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1. INTRODUCTION 

Magnetars are isolated neutron stars that have very high 
magnetic-field strengths inferred from the high spin down 
rate, typically greater than 10 14 G (e.g., Vasisht & Gotthelf 
1997; Kouveliotou et al. 1998; Olausen & Kaspi 2014, for 
a catalog of magnetars). Their X-ray luminosities are often 
greater than their spin-down power, and theorized, therefore, 
to be powered by the decay of intense internal magnetic 
fields (Duncan & Thompson 1992; Thompson & Duncan 1995, 
1996). The decay of the magnetic field may gradually build 
up stress in the crust, which can fracture it and/or twist the 
external magnetic field, causing short X-ray and soft-gamma- 
ray bursts (Feroci et al. 2001; Lenters et al. 2003; Gog lis 
et al. 2011) and sudden flux increases (Thompson et al. 2002; 
Beloborodov 2009). Historically, two classes of X-ray pulsars 
were thought to be magnetars: anomalous X-ray pulsars (AXPs) 
whose X-ray luminosities exceed the spin-down power, and soft 
gamma repeaters (SGRs) which show repeated soft gamma-ray 
bursts. However, distinction between the two classes has been 
significantly blurred (Gavriil et al. 2002; Kaspi et al. 2003). 
There are 26 magnetars discovered to date with various spectral 
and temporal properties (see Olausen & Kaspi 2014). 11 

Magnetar bursts have a variety of morphologies, including 
short (~ 100 ms) symmetric bursts, multiple peaked bursts, those 
with fast rises and longer decays, and some which exhibit very 


10 Lome Trottier Chair; Canada Research Chair. 

11 See the online magnetar catalog for a compilation of known magnetar 
properties: http://www.physics.mcgill.ca/~pulsar/magnetar/main.html. 


long extended “tails” (see, e.g., Woods et al. 2005; Gavriil et al. 
2006). Previous studies have suggested relationships between 
burst intensity and tail energetics (Lenters et al. 2003; Woods 
et al. 2005) and even the possibility of two distinct types of bursts 
(Woods et al. 2005). Burst spectra are generally described with 
thermal models (van der Horst et al. 2012) although there has 
been evidence for spectral features in some bursts (Gavriil et al. 
2002, 2011; Woods et al. 2005; Dib et al. 2009). 

The magnetar IE 1048.1—5937 is relatively active, often 
showing X-ray bursts and unstable timing behavior. Its spin 
period is 6.46 s, and the spin inferred surface magnetic-field 
strength is B = 4 x 10 14 G. In quiescence, it shows a 
spectrum which is well described with a blackbody plus power- 
law model having kT ~ 0.6 keV and T ~ 2.9 (Tam et al. 
2008). The distance to the source is estimated to be 9 kpc 
(Durant & van Kerkwijk 2006), which we use throughout 
this paper. Interestingly, it is the first AXP in which X-ray 
bursts were seen (Gavriil et al. 2002), and has shown several 
more bursts since then (Dib et al. 2009), hereby blurring the 
distinction between the AXP and SGR classes (see also Kaspi 
et al. 2003, for IE 2259+586). Another interesting property 
of IE 1048.1— 5937’s bursts is a possible emission feature at 
~ 1 3 ke V in the spectrum which was previously seen during its 
2002 burst in Rossi X-ray Timing Explorer (RXTE) data (Gavriil 
et al. 2002; Dib et al. 2009). Similar spectral features have been 
seen in X-ray bursts from other magnetars as well, all with RXTE 
(XTE J 1 8 1 0- 1 97, 4U 0142+6 1 ; Woods et al. 2005 ; Gavriil et al. 
2011 ). 

In this paper, we report on the spectral and temporal prop- 
erties of eight new bursts from IE 1048.1—5937 detected 
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serendipitously with the Nuclear Spectroscopic Telescope Array 
( NuSTAR ) in 2013 July. We describe the observations and data 
reduction in Section 2, and show the data analysis and results 
in Section 3. We then discuss the implications of the analysis 
results (Section 4) and conclude in Section 5. 

2. OBSERVATIONS 

NuSTAR operates in the 3-79 keV band, and is the most 
sensitive satellite in the ~ 10-79 keV band thanks to its unique 
focusing capability in this energy range. The energy resolution is 
400 eV at 10 keV (FWHM), and the temporal resolution is 2 /xs 
(see Harrison et al. 2013 for more details). Its excellent temporal 
and spectral resolutions together with the high sensitivity 
above ~10 keV are optimal for the study of burst properties 
of magnetars because such events are very brief (~ms) and 
sometimes spectrally very hard. 

IE 1048.1—5937 was observed with NuSTAR on 2013 
July 17-27 with a total net exposure of ~320 ks (Obs. 
ID 30001024001-7). A 70 ks joint XMM-Newton observation 
(Obs. ID: 0723330101) was conducted using the small window 
mode for MOS 1 /2 and the full frame mode for PN on 20 1 3 July 
22 to extend the spectral coverage down to ~0.5 keV where the 
thermal component is dominant. The source was not known to 
be in a particularly active state at the time of the observation. 
During the NuSTAR observations, eight X-ray bursts from the 
source were detected to our surprise, with one simultaneous 
detection in the XMM-Newton data. 

The NuSTAR data were processed with nupipeline 1.3.1 
along with CALDB version 20131007 using standard filters 
except for PSDCAL. We set PSDCAL=N0 in order to recover more 
exposure by slightly sacrificing the pointing accuracy. 12 We 
verified that the exposure increases and that the imaging, timing, 
and spectral analysis results are consistent with those obtained 
with PSDCAL=YES. The XMM-Newton data were processed with 
Science Analysis System 12.0.1 using the standard filtering 
process. We then further processed the event files for analysis 
as described below. 

3. DATA ANALYSIS AND RESULTS 
3.1. Burst Morphology 

In order to search for bursts, we extracted events in an aperture 
of 60" around the source position in the NuSTAR image, and 
produced barycenter-corrected light curves with a bin size of 
0.5 s. We then searched for time bins having a significantly larger 
number of counts compared with the persistent level which was 
extracted from the source region (and which is dominated by the 
persistent flux from the source) in 10 pulse periods (~64 s) prior 
to the time bin that is being searched. The average count rate of 
the persistent emission was 0.2 count s _1 in the 3-79 keV band 
within the R = 60" aperture (see Figure 1). We further verified 
that the count rate in the off-source region did not show any 
significant increase (e.g., due to a background flare) over any 
short time interval. We then calculated the Poisson probability 
of the observed count given the background rate for each time 
bin. To be considered significant, we required a >3er chance 
probability p after considering the number of trials. In total, we 
detected eight bursts with high significance, each one having 
p < 10~ I() . We also tried different bin sizes (e.g., 0.01-100 s), 
and found similar results, although the significance changes 


12 http://heasarc.gsfc.nasa.gov/docs/nustar/analysis/nustar_swguide.pdf 




Figure 1 . Top: 3—79 keV light curve for the NuSTAR observations with Modules 
A and B combined. Time intervals for pre-burst background extraction are shown 
in red hatched lines (see Sections 3.3 and 3.4). Insets show zoomed-in 3-79 keV 
NuSTAR and 0.5-10 keV XMM-Newton light curves (blue) around the bursts 
with burst numbers. Bottom: a 30 s light curve around burst 5 in the 3-79 keV 
band. T = 0 is set to the burst peak time, and red vertical lines show a one 
pulse period in which we characterize the burst morphology (see Equation (1) 
and Table 1). 

(A color version of this figure is available in the online journal.) 

depending on the binning, and some bursts are not significantly 
detected on some timescales. The bursts 2 and 7 are not detected 
on a long timescale (100 s), and we do not report their tail 
properties. 

We note that many time bins immediately after a burst 
were detected significantly above quiescence. The significantly 
detected time bins are likely those of the burst tail emission or 
the peaks of the source pulsations which are visible immediately 
after the burst (e.g., see Figure 1 bottom). However, there may 
be time bins with significantly higher count rates which are 
produced by an independent burst event. To investigate whether 
these significantly detected bins after a burst were independent 
events and not related to the burst tail or the pulse peaks, we 
proceeded as follows. We extracted a ~50 s light curve after 
a burst and characterized it with a decay model, specifically 
an exponentially decaying sine plus an exponential decay plus 
a constant. We then searched for time bins (AT = 0.05-1 s) 
having counts above the decay model with >3 ct confidence, 
and found none. Therefore, we conclude that there were only 
eight significant bursts during the observation. The NuSTAR 
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Summary of NuSTAR-detected Bursts 
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Notes. Uncertainties 

are at the 1 a confidence level, and upper limits are at the 90% confidence level. 



a Parameters for the short timescale li 

ght curve as 

defined in Equation (1). Tq is 

days since 

MJD 56490 (TDB). 




b Spin phase corresponding to To, where phase zero is defined at the pulse minimum (56490.3343345727 MJD), same as that for the 
timing analysis in Figures 2-4. 

c Time interval which includes 90% of the burst counts (the exponential functions in Equation (1)). 790 ’s for the rising and the falling 
function were calculated separately and then summed to obtain that for the burst. When only an upper limit is available for T r or Tf, we 
used the upper limit to calculate Tg o and show it without uncertainties. 

d Spectral parameters for a blackbody spectrum corresponding to Tgo'. blackbody temperature kT and bolometric luminosity Lgo- 


light curve is shown in Figure 1, and the burst properties are 
listed in Table 1. The burst light curves all exhibit rises and 
decays of a few seconds with relatively long tails (~ks). 

After we identified the bursts, we extracted events in a 
one period interval that contains the burst in order to further 
characterize the temporal properties of the bursts on very short 
timescales. Note that we included one more period for the bursts 
that occurred late in pulse phase so as not to miss the falling 
tail. We fit the time series to an exponentially rising and falling 
function 


F(f) = 


Ae «-To)/T, + Cx 

(A - C 2 )e- (f - 7b)/7 ’ f + Ci+C 2 


t < Tq, 
t ^ To, 


( 1 ) 


where A is the amplitude, Tq is the burst peak time, T r is the 
rising time, Tf is the falling time, and C' 1.2 are constants. We 
note that the decay of a burst has typically been modeled with 
a double exponential function (e.g., Gavriil et al. 2011), but 
here we replace the second exponential with a constant (C 2 ) 
because this suffices for describing our data in the chosen 
time span, which is much smaller than the decay constant of 
any second exponential. Since the timescales are very short, 
and there are only a few events in each ~ms time bin, we 
used maximum-likelihood optimization. Furthermore, we used 
events in the whole detector because having a well sampled time 
series is important in the fitting. We modeled the background 
with another constant, C\ . The results of the fitting are presented 
in Table 1 . 

We note that the observed count rate is smaller than the 
incident rate for the brightest burst due to detector deadtime 
(Harrison et al. 2013). Since the maximum count rate for 
burst 4 is comparable to the maximum count rate that the 
NuSTAR detectors can process ("400 counts s~ 1 per module), 
we consider the effect of deadtime in order to calculate the 
incident count rate, R[, via the following relation: 


1 - R 0 r D 

Here R a is the observed count rate, and r n is the detector dead 
time (~2.5 ms) for each observed event. The incident peak count 


rates are higher, and the rising and the falling times are smaller 
than the observed values in Table 1, but within a factor of "2 of 
the true values. For example, the maximum observed count rate 
for burst 4 is ~200 counts s -1 per module, and the incident rate 
is estimated to be ~400 counts s' 1 per module for this burst. 


3.2. Timing Analysis 

In order to measure the pulse period, we extracted events 
within radii of 60" and 32" for NuSTAR and XMM-Newton, 
respectively, and applied a barycenter correction to the events. 
We then subdivided the total observations into 5 1 and 48 sub- 
intervals so that each sub-interval has ~ 1 200 and ~2400 events 
to have at least ~20 counts in a phase bin (Acp — 1/32) 
for NuSTAR (3-79 keV) and XMM-Newton (0.5-10 keV), 
respectively. The typical duration of each sub-interval is ~15 ks 
and ~1 ks for NuSTAR and XMM-Newton, respectively, but 
varies depending on the source luminosity. We then fit the pulse 
profile to a Gaussian plus constant function to measure the phase 
of each profile. The Gaussian plus constant function describes 
the pulse profiles well, and we show the measured phases in 
Figure 2 and examples of pulse profiles in Figures 3 and 4. 
We note that cross-correlating the pulse profiles gives similar 
results. The phases are fitted to a linear function cp(t) — 4>o+ ft, 
where <po is the reference phase, and /is the frequency. We did 
not include the frequency derivative because it is not required. 
In fitting, we ignored ~ 10-90 ks of data after the bursts because, 
interestingly, there is a relatively large phase shift immediately 
post-burst which contaminates the result (see Figure 2). From 
the analysis, we found the period to be 6.46168155(6) s, but 
were not able to constrain the period derivative well. 

In order to see if the large phase shift during the bursts (see 
Figure 2) is related to the rotation of the star, we measured the 
shift in different energy bands ( NuSTAR and XMM-Newton) and 
found that it is much smaller in the soft band than in the hard 
band (e.g., see T ~ 6 days in Figure 2). We also verified that 
the energy dependence of the post-burst phase shifts is observed 
in the NuSTAR and the XMM-Newton data individually, and the 
large phase shift at ~6 days in the NuSTAR data was seen in the 
XMM-Newton hard-band data (4-10 keV) as well. These imply 
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Figure 2. Measured timing residuals for P = 6.46168155 s for the NuSTAR 
data in the 3-79 keV band (black) and the XMM-Newton data in the 0.5-10 keV 
band (red). The best-fit function is shown in a blue horizontal dashed line and 
the burst times are shown in green vertical dotted lines. 

(A color version of this figure is available in the online journal.) 

that the phase shift is not due to the rotation of the source; if it 
were, we would expect the shift to be independent of energy. 

Figures 3 and 4 present pre-and post-burst pulse profiles for 
IE 1048.1—5937 from NuSTAR and XMM data, respectively. 
The pre- and post-burst pulse profiles are qualitatively similar 
except for the phase shifts discussed above. Some “spiky” 
features appeared in some post-burst profiles (e.g.. Figure 3), 
which are likely to be statistical fluctuations. However, we note 
that another hard peak at phase ~0.3 in the 3-10 keV band 
seemed to appear in the XMM-Newton pulse profile immediately 
post-burst (see Figure 4(b)). 


3.3. Burst Spectroscopy 

Next, we focus on the spectra of the bursts and their tails. 
The spectrum of the full data will be presented elsewhere (R. F. 
Archibald 2014, in preparation). In order to characterize a burst 
spectrum, we calculated the Tgo for each event; these are shown 
in Table 1 . We note that deadtime is likely to have affected the 
burst light curves (Section 3.1) since we used unbinned events 
and did not correct for the deadtime for individual events. A 
spectrum, integrated over a time interval, is less affected since 
the deadtime effect is corrected for every 1 s time bin, which is 
precise enough unless the spectral shape rapidly changes within 
the 1 s time bin. We assume that the spectral shape (i.e., kT) of 
the source did not change significantly in one second, and hence 
that the deadtime effect is properly corrected. 

We extracted source events in a circle with radius 60" in the 
time interval T 90 for each burst. Backgrounds were extracted 
from the pre-burst time intervals shown in Figure 1 (red hatched 
regions). The pre-burst spectra include photons up to ~ 15 keV 
above the background and are well described with a power law 
plus blackbody or two blackbody models, both having 3-79 keV 

flux of 4 x 1 0 12 erg s' 1 cm 2 . The luminosity for the assumed 

two blackbody model is 1 x 10 35 erg s' 1 (see R. F. Archibald 
2014, in preparation for more detail). Note that the extraction 
times are very short, so background contamination is very small 
(~0.2 cps including the persistent source emission). 

We fit the spectrum to an absorbed blackbody and an absorbed 
power law using the C statistic (cstat in XSPEC 12.8.1) in the 
3— /-max keV band, where £ nlax was between 30 and 50 keV 
depending on the source flux. We verified that the fit results 
did not change over a broad range of the upper energy limit. 
There were not many events in each spectrum, and we found 
that both spectral models provided good fits. Since NuSTAR is 
not sensitive in the energy band below 3 keV, we were not able to 



Figure 3. Immediate pre- (black) and post-burst (red) NuSTAR pulse profiles for bursts 1, 3, 4, 5, 6, and 8 in the 3-79 keV band. Tgo burst intervals were not included 
in making these profiles. Backgrounds were subtracted, and the integration times were chosen so that each pulse profile had ~1200 counts. Burst numbers are shown 
in the plots for reference. Vertical lines in the plots are the phases corresponding to the 7o in Equation (1) for each burst. Note that the pre-burst profiles of burst 5 and 
6 include significant tail emission of burst 4 and 5, respectively. 

(A color version of this figure is available in the online journal.) 
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Figure 4. XMM-Newton pulse profiles of pre- and post-burst time intervals of burst 4 in the 0.5-3 keV (black) and 3-10 keV (red) bands for the pre-burst (a), 
immediately post-burst in a time interval of 2-2000 s after the burst peak time measured with NuSTAR (b), and late post-burst (c) time intervals. Backgrounds were 
subtracted, and the integration times were 2 ks. The hard-band profiles were normalized to the soft-band ones for each plot. 

(A color version of this figure is available in the online journal.) 


measure /V H , and so we froze it at the value measured previously 
(0.97 x 10 22 cm' 2 ; Tam et al. 2008). 

We show the best-fit parameters for the blackbody model in 
Table 1 . Note that the numbers of events in Tgo for the bursts 2, 
3, 4, and 7 were smaller than ~10, and the spectral parameters 
were not reasonably constrained. We also tried to fit the data 
using the usual x 2 method with Churazov weighting (Churazov 
et al. 1996) after grouping the spectra to have at least 5 counts 
per bin found that the results are consistent with those obtained 
using the C statistic. 

Although single component models provide good fits, we 
also tried to fit the spectra to double blackbody models, 
since magnetars sometimes show double blackbody spectra 
during bursts (e.g., Lin et al. 2012). The second blackbody 
component is not statistically required for any bursts in our 
data. However, it is possible that there was an undetected high- 
temperature blackbody component as was seen by, e.g., Lin 
et al. (2012). Therefore, we measure a luminosity upper limit 
for a high-temperature blackbody component. We first froze the 
parameters of the single blackbody fit at the values obtained 
above, added another high-temperature blackbody, and scanned 
the parameter space of the high-temperature blackbody using 
the steppar command in XSPEC. The luminosity upper limit 
increased with kT as expected since NuSTAR becomes less 
sensitive as kT goes up. We measured the 90% luminosity upper 
limit to be 0.2-8 x 10 38 erg s' 1 for kT = 15 keV, for example. 

We also set an upper limit on flux of a possible lower 
temperature blackbody, similar to that seen in the bursts of 
SGR 1900+14 (Israel et al. 2008), where the authors found 
a low-temperature blackbody component with kT as low as 
~2 keV having comparable luminosity with the high temper- 
ature component. We followed the same procedure that we 
did above for the high-temperature component but with kT = 
2 keV, and found that the luminosity upper limits are 2-10 x 
10 36 erg s' 1 , always an order of magnitude smaller than the burst 
luminosities in Table 1 . 

3.4. Tail Spectroscopy 

We also characterized the post-burst tail emission for the 
bursts with a significant tail (bursts 1, 3, 4, 5, 6, and 8). In order 
to minimize the effect of the burst, we removed the Tgo intervals 
in this analysis. We extracted events from a circular region 
of a radius 60" centered on the source, and the background 
from a source-free region on the same detector in a radius 90" 
in the NuSTAR data. Since the source spectrum may change 
significantly during a burst, we subdivided each burst into 


several sub-intervals such that each interval has ~200 counts, 
in order to have good time resolution and to be statistically 
sensitive to the change in the spectral parameters in time. The 
integration times were very short at early times after a burst, 
and we could not collect enough background events. Therefore, 
for those sub-intervals we used a longer background exposure 
(1 ks). 

We fit the spectrum for each sub-interval to an absorbed black- 
body and an absorbed power-law model using the C statistic in 
XSPEC. The fit ranges were 3-£ max , where £ max was 20-50 keV. 
Both blackbody and power-law models are acceptable, although 
blackbody models provide better fits in general, giving fit statis- 
tics smaller by ~10 on average (~200 dof). We show the results 
for the blackbody fit in Figures 5(a)-(c). Note that the persistent 
emission is included in the spectra in this case. 

We then tried to characterize the spectral evolution after 
removing the persistent emission. The persistent level was 
extracted from the source region but in the pre-burst time 
intervals shown in Figure 1. The source seemed to return to 
the persistent level ~2-3 ks after a burst, and thus we analyzed 
only the first ~2 ks after a burst. 

We find that the spectral parameters and their evolution are 
similar to those of the spectra with the persistent emission. This 
is expected because the burst tails dominate over the persistent 
emission during the time intervals. 

The bursts have very long tails (~ ks) as is often seen in 
magnetar bursts (e.g., Gavriil et al. 2011). We have measured 
the timescales of the decaying tails by fitting the light curve. 
However, as we show above (Figure 5), the spectrum evolves 
significantly over a relatively short timescale, and hence the 
measured light curves obtained using energy-integrated count 
rates will be different for different telescopes because of 
differences in the energy responses. In order to estimate the 
decay timescales in an instrument-independent way, we directly 
measure the spectral evolution of the bursts. We fit the time 
evolution of the spectral parameters for the tail spectra including 
the persistent emission with a power-law decay, 

S = S 0 r a + Si, (2) 

where S is the spectral parameter, So is the value of the parameter 
at t = 1 s, a is the decay index, and S) is a positive constant 
corresponding to the persistent emission. We show the results in 
Table 2. Note that we do not show the results for burst 3 because 
we were not able to constrain its decay with the given data. 

We note that there might be some less significant bursts in the 
tails. For example, bursts 2 and 7 occurred at ~600 s and ~200 s 
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Figure 5. Time evolutions of the spectral parameters of the blackbody model for the persistent plus tail spectra and kT vs. Lbb for tail spectra of the six bursts (bursts 
1, 3, 4, 5, 6, and 8). (a) Blackbody temperature, kT, (b) bolometric luminosity, Lbb^ (c) blackbody area, and (d) kT vs. Lbb- 
(A color version of this figure is available in the online journal.) 


Table 2 

Parameters for the Spectral Evolution of the Tails 


Burst 

So 

(10 35 erg s~*) 

<*L 

a kT 

Si 

(10 35 erg s” 1 ) 

1 

290(60) 

0.81(4) 

0.28(1) 

0.69(7) 

4 

90(50) 

0.8(1) 

0.57(5) 

0.84(6) 

5 

640(90) 

0.82(3) 

0.27(1) 

0.5(1) 

6 

560(120) 

0.87(5) 

0.29(1) 

0.7(2) 

8 

420(150) 

1.03(9) 

0.28(4) 

1.21(8) 


Note. Uncertainties are at the la confidence level. 


into the tail of burst 1 and 6, respectively. Also, there seems to 
be an increase in kT and Lbb at ~600 s after burst 5 although we 
did not find any significant burst. Undetected or less significant 
bursts may bias a to a smaller value. Also note that we assumed 
that the luminosity at late times can freely vary in the fit, and 
the best-fit values for Si are 0.5-1. 2 x 10 35 erg s” 1 , sometimes 
smaller than the pre-burst source luminosity. If we freeze Si at 
the pre-burst value of 1 x 10 35 erg s 1 , a becomes 0. 9-1.0. 

We also analyzed the XMM-Newton data to characterize 
precisely the soft-band (0.5-10 keV) spectrum of the burst 4 
tail which is the only burst detected simultaneously with XMM- 
Newton and NuSTAR. We extracted source and background 
events in circles with radius 32" for an exposure of ~400 s 
excluding the Tgo interval, where the background region was 
~200" north of the source. The XMM-Newton count rates 
were 1.2/1. 3 and 4.5 cps for MOS1/2 and PN, respectively. 
Such count rates will result in spectral distortion of <1% and 
~2.5% for MOS1/2 and PN, much smaller than the statistical 
uncertainties we obtain below. Hence, we used all the MOS 1 /2 
and PN data. We grouped the spectrum to have at least 20 events 
per spectral bin. We then extracted NuSTAR events for the same 


time interval, and binned the spectrum to have 20 events per bin, 
and jointly fit the XMM-Newton and the NuSTAR data, holding 
Ah fixed at 0.97 x 10 22 cm” 2 . Note that the persistent emission 
was not removed in these spectra. 

Single component models were not acceptable with / 2 / 
dof of 583/102 and 122/102 for a blackbody and a power 
law, respectively. A double blackbody model was acceptable 
with x 2 /dof = 127/100. But a blackbody plus power-law 
model provided a better fit (/ 2 /dof = 98/100), and adding 
one more blackbody slightly improved the fit (x 2 /dof = 
96/98) although it was required only with low significance 
(/•’-test probability 40%). Nevertheless, the parameters for the 
soft spectral component (k 7] = 0.6±0.2keV andT = 2.5±0.4) 
of the power law plus two blackbody model were similar to those 
of the quiescent spectrum of the source (Tam et al. 2008), and 
the hard component ( kT h — 3.2 ± 0.4 keV) is similar to that of 
the spectrum with the persistent emission removed (see below). 
When we let Ah vary, we obtain similar results as the above 
with A h = 1.2(4) x 10 22 cm” 2 . 

We also tried to characterize the persistent-emission-removed 
spectrum of the combined data. Here, the persistent level was 
extracted from the pre-burst time intervals in the source region 
for both NuSTAR and XMM-Newton. A blackbody model was 
able to fit the data well (kT = 3.0±0.2 keV, x 2 /dof= 96/104) 
but a power-law model was not as good (x 2 /dof = 128/104). 

3.5. Spectral Feature 

We searched for the spectral feature at ~13 keV previously 
observed from this source’s burst spectra (Gavriil et al. 2002, 
2011; Dib et al. 2009). As the feature was observed in the <2 s 
spectra in the past, we checked if a line feature was statistically 
required in the burst spectra (see Section 3.3) but found that it is 
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Table 3 

Best-fit Spectral Parameters for the Spectrum of the First 150 s after To for Burst 5 


Model 3 

Ah b 

(10 22 cnU 2 ) 

kT/r 

(keV) 

Lbb/F pl c 

E G d 

(keV) 

<^G 

(keV) 

Nq 

(photon cm -2 s -1 ) 

X 2 /dof 

BB 

0.97 

4.1(2) 

30(2) 




70/42 

BB + Gauss 

0.97 

3.8(2) 

25(2) 

13.2(3) 

1.3(3) 

2.1(5) x 10^ 3 

39/39 

PL 

0.97 

0.83(7) 

12(1) 




108/42 

PL + Gauss 

0.97 

1.0(1) 

0.8(1) 

12.7(2) 

1.7(3) 

3.4(6) x 10^ 3 

44/39 


Notes. Uncertainties are at the ltr confidence level. 
a BB: blackbody, PL: power law. Gauss: Gaussian line profile. 
b Frozen. 

c Bolometric luminosity in units of 10 35 erg s 1 for the blackbody model, and the absorption corrected 3-79 keV 
flux in units of 1 0 3 erg s 'em 2 for the power-law model. 

d Subscript G is used for the Gaussian parameters: line energy (Eq), line width (<tgX and normalization (Na)- 



Figure 6. Blackbody plus Gaussian line fit for the 150 s NuSTAR spectrum of 
burst 5. 

(A color version of this figure is available in the online journal. 1 


not. However, as we increased the integration time, we started 
seeing an enhancement in counts at ~13 keV for burst 5. 

In order to see if a line feature is statistically required 
for the bursts and their tail emission, we first extracted a 
spectrum from the first 150 s of the brightest burst (burst 5). 
The background was extracted from the pre-burst interval in 
the source region. We then fit the spectrum in the 3-30 keV 
band with continuum models (blackbody, power law, and 
bremsstrahlung), and found that none could describe the data 
well with / 2 /dof of 70/42 (blackbody), 108/42 (power law), 
and 147/42 (bremsstrahlung). We also tried combinations of 
the continuum models, but none improve the fit significantly 
with x 2 /dof being 67/40, 108/40, 67/40, 67/40 for a double 
blackbody, a double power-law, a blackbody plus power-law, 
and a blackbody plus bremsstrahlung models, respectively. 
In particular, in all these two-component trials, the best-fit 
parameters of the second component are trivial; the amplitude 
is not constrained at all, and kT or T become ridiculously small 
or large. However, adding a Gaussian line to a blackbody model 
provided an acceptable fit (/ 2 /dof = 38/40). We show the 
spectrum and the blackbody plus Gaussian fit in Figure 6, and the 
results are presented in Table 3. We repeated the same procedure 
for spectra in different time intervals (e.g., 100 s, 200 s, and 
300 s), and found that the fit improves significantly (A/ 2 > 20) 
by adding a Gaussian line to the blackbody model, from which 


we draw the same conclusion as we did with the 150 s spectrum. 
We conducted the same analysis with the other bursts and found 
that a blackbody model alone was able to describe their spectra 
well (for example, / 2 /dof = 39/35, p ~ 0.1 for burst 6). 

We used simulations in order to calculate the true significance 
of the feature in burst 5 and its tail. The spectrum of the source 
following a burst evolved on a timescale of 150 s. Therefore, 
we conducted simulations with temporally evolving blackbody 
spectra. We first divided the 150 s spectrum of burst 5 into 
five time intervals as shown in Figure 5 and the Tijo interval, 
and fit the six spectra to blackbody models after removing the 
energy range of the feature (11.5-14 keV). We then simulated 
the six spectra using fake it in XSPEC with the response and 
background for the actual data, and merged them into one 
spectrum that represents a 150 s spectrum. For each simulation, 
we first fit the spectrum to a blackbody, and then calculated 
the improvement of the fit (Ax 2 ) by adding a Gaussian line. 
For the latter fit, we scanned the energy range of the spectrum 
(4-22 keV) with a step of 1 keV for the initial value of the line 
energy because the narrow feature can make the fit fall in a local 
minimum around the initial value. We counted the number of 
occurrences in which the improvement of the fit in a simulation 
was larger than that seen in the actual data fit (Ax 2 = 31), 
and found that this did not occur in 10,000 simulations. The 
significance of the feature would be even higher if we scan only 
a smaller energy range (e.g., 11-15 keV) for the Gaussian line 
energy guided by the previous RXTE measurements. 

We found that a blackbody model fits the spectra of the 
feature-removed 150 s data as well as the above simulations. 
Therefore, we also conducted simulations with a single black- 
body continuum model, and found that improvement of fit mea- 
sured by Ax 2 by adding a Gaussian line model was always 
smaller than 31 in 10,000 simulations. 

We investigated if the feature exists in the spectra without the 
burst. We removed the Tgo interval from burst 5 and extracted 
tail spectra for 100 s, 150 s, 200 s, and 300 s exposures. We 
then applied the same fitting procedure and found that the 
same conclusion is valid; for example, single or combinations 
of continuum models do not fit the 150 s spectrum well 
with x 2 /dof of 75/39 (blackbody), and 164/39 (power law) 
while adding a Gaussian line to a blackbody improves the 
fit significantly, making it acceptable (x 2 /dof = 42/36). The 
best-fit Gaussian line parameters are statistically consistent with 
those we obtained above. 

We checked if the spectral feature shifts in energy with phase 
as seen in the magnetar SGR 0418+5729 (Tiengo et al. 2013). 
We produced a two-dimensional energy versus phase image as 
shown in Tiengo et al. (2013), and found no evidence of a shift. 


7 



The Astrophysical Journal, 790:60 (9pp), 2014 July 20 

This may be due to the paucity of counts in our case. We also 
tried to see if the feature is more prominent in some phase 
intervals than the others, and could not draw any statistically 
significant conclusion. 

4. DISCUSSION 

We have reported the first detection of X-ray bursts from 
magnetar IE 1048.1—5937 with a focusing hard X-ray tele- 
scope. Table 1 summarizes the temporal and spectral properties 
of the eight bursts detected with NuSTAR. The short timescale 
behaviors of the bursts are different from burst to burst, having 
T t and 7f of milliseconds to seconds. All the bursts have rel- 
atively long tails (see Figure 5) and their spectra evolve with 
time. We also found that their peak times are random in phase. 

4.1. Temporal Properties of the Bursts and their Tails 

Magnetar bursts, including those we have observed for 
IE 1048. 1—5937, generally exhibit a fast rise and fall, which has 
been suggested to result from crustal fracture (Thompson et al. 
2002) and/or magnetic reconnection (Lyutikov 2003). The pulse 
phase of the burst peak was previously used to reconstruct the 
event location and constrain the emission geometry (Woods et al. 
2005). Interestingly, the bursts appear to occur at random pulse 
phases (Table 1), and bursts 2 and 5 may be emitted far from 
the magnetic axis of the star. Note, however, that the location 
reconstruction assumes emission beaming along the magnetic 
field lines, which may not be realistic. A quasi-isotropic burst 
could be visible from any direction unless it is eclipsed by the 
neutron star. The probability of eclipse is significantly reduced 
by gravitational light bending, which makes »3/4 of the star 
visible to a distant observer (Beloborodov 2002). 

A burst can have tail emission which is produced by the 
residual heat of the crustal fracturing (Lyubarsky et al. 2002) 
or bombardment of the stellar surface by the magneto spheric 
particles (Beloborodov 2009). The tail emission lasts much 
longer than a burst, and allows us to sample many full rotations. 
Studying pulse profiles during burst tails may provide insight 
into the tail emission region. Interestingly, we find clear phase 
shifts in some post-burst profiles (e.g., bursts 4 and 5, see 
Figures 2 and 3), which are not caused by rotation (Section 3.2). 
We conclude that the source of tail emission, e.g., the hot spot 
produced by the burst, is slightly shifted in longitude relative to 
the source of persistent emission. 

The pulse profiles measured by XMM-Newton (Figure 4) also 
hint at the appearance of a new hot spot. The immediate post- 
burst pulse profile shows an additional peak at phase ~0.3 
(Figure 4(b)), which may correspond to the new hot spot. Its 
emission added to the pulsed persistent source results in a phase 
shift, which returns to zero as the luminosity of the new hot spot 
decays and the persistent source dominates again (Figure 2). 

4.2. Burst Spectra and Their Evolution 

Interesting correlations between burst spectral properties have 
been seen in some magnetars. For example, Gdgiis et al. (2001) 
reported an anti-correlation between fluence and hardness ratio 
in the bursts of SGR 1806-20 and SGR 1900+14; however, 
Gavriil et al. (2004) observed the opposite trend between the 
same properties in the bursts of IE 2259+589. While detailed 
studies of such correlations here are not possible due to the 
limited number of bursts and small statistics in each burst, 
we note that there seems to be a hint of a positive correlation 
between L 90 and kT (see Table 1). 


An et al. 

The source flux returns to the pre-burst value on the kilo- 
second timescale. Such a short timescale suggests that the 
crustal heating was not deep, as the large specific heat of the 
deep crust would hold the energy longer and delay the decay 
(Kouveliotou et al. 2003). The spectrum of tail emission at 
~ 10-1000 s can be fitted by a blackbody or a power-law 
model; the blackbody model provides a better fit when the 
persistent emission is subtracted. The evolution of the blackbody 
luminosity Lbb can be described as a power law with the decay 
index of 0.8-0. 9 (Table 2 ). 13 Similar flux decay indices were 
measured in the 2004 bursts of IE 1048—5937 (Gavriil et al. 
2006, 0.82 ± 0.05) and XTE J1810-197 (Woods et al. 2005, 
0.92 ± 0.02). The indices are not far from the crustal cooling 
model with Lbb oc t~ 2 ^ 3 (Lyubarsky et al. 2002), which was 
previously used for later time evolution ( t > 10,000 s). The 
/-bundle untwisting (Beloborodov 2009) is usually invoked on 
even longer timescales, which are associated with the resistive 
evolution of the magnetosphere. 

Note the increase in blackbody area (Figure 5(c)) following a 
burst. We argue that this is an observational effect rather than a 
real increase in the physical size. If a burst produced a local hot 
spot as we argued above, the measured size of the blackbody 
area would represent the small burst spot at early times, and the 
large persistent spot later. In order to disentangle the effects of 
the hot spot and the persistent spot, we also measured the size 
evolution with the persistent-emission-removed spectrum. The 
tail was detected significantly above background only for the 
first ~ 1 ks, and we found that the blackbody area evolved in a 
similar manner to the trend of the first ~1 ks in Figure 5; no 
clear increase in the size was observed. 

In magnetar cooling scenarios, an excited magnetar cools 
by emitting thermal photons and/or untwisting a /-bundle. In 
both cases, a correlation between flux and hardness is expected 
(Thompson et al. 2002; Lyutikov 2003; Ozel & Giiver 2007). 
Such a correlation has been generally observed in magnetars’ 
short-term and long-term cooling (e.g., Woods et al. 2005; 
Gavriil et al. 2006; Scholz et al. 2011; An et al. 2013) with 
some exceptions (e.g., An et al. 2012; Kaspi et al. 2014). 
We investigated if such a correlation exists in the tail spectra 
and found a clear correlation between kT and log 10 LBB (see 
Figure 5(d)), implying that the burst tails also exhibited a 
correlation between flux and hardness. Whether the origin of 
this correlation is the same as for that of the long-term cooling 
is not clear. 

In principle, tail radiation could be the burst “echo” produced 
by dust scattering around the magnetar (Tiengo et al. 2010). In 
this scenario, the tail spectrum should soften with time (a result 
of the scattering cross section being smaller at higher energies); 
whether or not this is consistent with the observed spectral 
evolution is unclear. Furthermore, the observed increases in the 
pulsed flux immediately after the bursts cannot be produced 
by dust scattering, suggesting that the tail is emitted by the 
magnetar itself. 

4.3. Spectral Feature 

We find a spectral feature at ~13 keV in the tail emission. 
Spectral features at a similar energy were previously observed 
from the burst and tail spectra of several sources, but only 
with RXTE (IE 1048.1-5937, XTE J1810-197, 4U 0142+61; 


13 The true decay may be somewhat steeper when the persistent emission of 
10 35 erg s _1 is subtracted; besides the measured decay index might be affected 
by undetected low-significance bursts during the tail emission. 
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Gavriil et al. 2002, 2006, 2011; Woods et al. 2005; Dib et al. 
2009). The line energy we found is similar to those previously 
reported. We note that this is the first detection of the feature with 
an instrument other than RXTE, which demonstrates the effect is 
not instrumental. We note that the line flux we measured is only 
~10% of that previously reported (Gavriil et al. 2002), which 
could be simply due to the long integration time we used (150 s 
versus 1 s). 

If we interpret the line-like emission as an electron cyclotron 
feature, the line energy implies that the magnetic field strength 
is ~1() 12 G in the emission region. A magnetic field strength 
of ~2 x 10 15 G could be inferred if we interpret the feature 
due to a proton cyclotron emission, and the line width implies 
a ~ 1 0% change in the magnetic field strength in the emission 
region. The three detections of the feature from IE 1048.1—5937 
(Gavriil et al. 2002; Dib et al. 2009, and this work) have similar 
properties (e.g., line energy and width), and probably originated 
from the same region of the star. If so, it is unlikely that a 
physical structure can be sustaining at a height of ~70 km 
from the stellar surface (where B ~ 10 12 G) to power bursts and 
line features multiple times over a decade. However, a strong 
multipolar magnetic field (10 15 G) near the surface in a volume 
of ~1 km 3 has enough energy to power multiple bursts and 
the line feature, suggesting that the line feature could be from 
proton cyclotron emission. Nevertheless, an interesting question 
is how different sources which may have different magnetic 
fields (IE 1048.1-5937, XTE J1810— 197, and 4U 0142+61) 
with different spectral and temporal properties show similar 
features at such similar energies. 

5. CONCLUSIONS 

We presented detailed spectral and temporal analyses of 
eight X-ray bursts from magnetar IE 1048.1—5937 detected 
with NuSTAR, one of which was simultaneously observed with 
XMM-Newton. The bursts exhibited a fast rise and decay with 
Tgo intervals of 1-4 s, and their spectra can be described with 
single blackbody models having kT ~ 6-8 keV for the Tgg 
intervals. All the bursts showed tail emission which can be 
described with temporally relaxing blackbody models. The flux 
relaxations of the tail emission followed a power-law decay 
having decay indices 0. 8-0.9. We confirm the existence of an 
emission feature at ~13 keV observed in burst and tail spectra 
of the source in the past. Finally, we note that similar spectral 
features at a similar energy have been seen in bursts of several 
magnetars with different spectral and temporal properties. This 
requires further theoretical interpretations. 
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